source('00_util_scripts/mod_seurat.R')

key_cytokine <- c('Ccl20','Tslp','Flt3l','Csf2','Tnf','Il6')

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg')

# import ---------
mex.path <- list.files('mission/FPP/zhang21sa', 'gz', full.names = T)

sobj <- mex.path |> read_geo_supp('(PBS|SA)_(HFD|SD)')

sobj <- sobj |>
  PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj |> VlnPlot('mito.ratio')

sobj <- sobj |>
  filter(mito.ratio < 20)

sobj %<>% quick_process_seurat()

micerna <- celldex::MouseRNAseqData()

sobj %<>%
  mark_cell_type_singler(micerna, new_label = 'mmur.main')

sobj |> DimPlot(group.by = 'mmur.main', cols = DiscretePalette(36))

DotPlot(sobj, c('Krt14','Krt15','Krt1','Krt10'))

sobj %<>% mutate(mmur.main = ifelse(seurat_clusters == 9,
                                  'Keratinocytes', mmur.main))

sobj |> DimPlot(group.by = 'mmur.main', cols = DiscretePalette(36))

sobj |> write_rds('mission/FPP/zhang21sa/zhang21sa.rds')

sobj <- read_rds('mission/FPP/zhang21sa/zhang21sa.rds')

sobj |>
  DotPlot(kegg_mva, cols = 'RdYlBu')

# KC --------
sobj.kc <- sobj |>
  filter(mmur.main == 'Keratinocytes')

sobj.kc |> DotPlot(kegg_mva, group.by = 'orig.ident', cols = 'RdYlBu') +
  RotatedAxis()

sobj.kc %<>%
  quick_process_seurat(skip_norm = T)

sobj.kc |>
  filter(orig.ident == 'PBS-SD') |>
  DotPlot(c(kegg_mva, 'Trpv3', late.kc), cols = 'RdYlBu') +
  RotatedAxis() +
  labs(title = 'Normal diet PBS mice skin KC: MVA pathway & Trpv3 expression',
       x = 'gene', y = 'cluster')

sobj.kc |>
  filter(orig.ident == 'SA-SD') |>
  DotPlot(c(kegg_mva, 'Trpv3', late.kc), cols = 'RdYlBu') +
  RotatedAxis() +
  labs(title = 'Normal diet SA mice skin KC: MVA pathway & Trpv3 expression',
       x = 'gene', y = 'cluster')

leiden.list <- sobj.kc$seurat_clusters |> unique()

savpbs.kc <- leiden.list |>
  map(\(x)FindMarkers(sobj.kc, group.by = 'orig.ident', subset.ident = x,
                      ident.1 = 'SA-SD', ident.2 = 'PBS-SD') |>
        mutate(cluster = x) |> as_tibble(rownames = 'gene'),
      .progress = T) |>
  list_rbind()

savpbs.kc |>
  filter(gene %in% kegg_mva) |>
  ggplot(aes(gene, cluster, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  labs(title = 'MVA pathway in KC: SA vs PBS',
       subtitle = 'GSE150729 (48 hpi)')

# T cell ---------
sobj.tc <- sobj |>
  filter(mmur.main == 'T cells')

sobj.tc %<>% quick_process_seurat(skip_norm = T)

sobj.tc |>
  DotPlot(map(seurat_markers, str_to_title))

sobj.tc |>
  DotPlot(c(t.submarker),
          cols = 'RdYlBu') +
  RotatedAxis()

immgen <- celldex::ImmGenData()

sobj.tc |>
  mark_cell_type_singler(immgen, fine_label = T)

sobj.tc <- sobj.tc |>
  mutate(type.fine = case_when(seurat_clusters == 5 ~ 'Th2 cells',
                               seurat_clusters == 1 ~ 'CD8 T cells',
                               seurat_clusters == 2 ~ 'Treg cells',
                               .default = 'gd Th17 cells'))

sobj.tc |> DimPlot(group.by = 'type.fine') +
  labs(title = 'T cell subsets in mouse skin',
       subtitle = 'GSE150729')

t.frac.2dpi <- sobj.tc |>
  as_tibble() |>
  filter(str_detect(orig.ident, 'SD')) |>
  calc_frac_conf_on_grouped_count(orig.ident, type.fine)

t.pval.2dpi <- sobj.tc |>
  as_tibble() |>
  filter(str_detect(orig.ident, 'SD')) |>
  test_on_grouped_count(orig.ident, type.fine)

t.frac.2dpi |>
  ggplot(aes(type.fine, fraction, color = orig.ident)) +
  geom_col(position = 'dodge2', fill = 'white') +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(width = 0.9), width = .5) +
  theme_pubr() +
  scale_color_hue(direction = -1) +
  labs(title = 'T cell subsets fraction change in SA skin infection',
       subtitle = 'GSE150729 (48 hpi)', x = 'cell type')

t.frac.2dpi |>
  left_join(t.pval.2dpi, join_by(type.fine == subtype)) |>
  write_csv('4dpi.T.subset.fraction.csv')

## T marker expression --------
t.mark.cura <- read_delim('t.subset.mark.curated.txt')

t.sublist <- sobj.tc$type.fine |> unique()

Idents(sobj.tc) <- 'type.fine'

t.mark.fc <- t.sublist |>
  map(\(x)sobj.tc |> FindMarkers(group.by = 'orig.ident',
                                 ident.1 = 'SA-SD', ident.2 = 'PBS-SD',
                                 subset.ident = x,
                                 logfc.threshold = 0) |>
        as_tibble(rownames = 'gene') |>
        mutate(cluster = x), .progress = T) |>
  list_rbind() 

t.mark.fc |>
  inner_join(t.mark.cura) |>
  filter(!str_detect(cluster, 'reg'), !str_detect(cluster, 'h1 ')) |>
  ggplot(aes(gene, cluster, color = avg_log2FC, size = -log(p_val_adj))) +
  geom_point() +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,3)) +
  facet_wrap(~cluster, scales = 'free', ncol = 1) +
  labs(title = 'T cell subset marker change: 4dpi SA vs PBS') +
  theme_jpub +
  rotate_x_text(45)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/R3.Q10.figI.T.subset.marker.4dpi.EF.vs.PBS.csv')

publish_pdf('mission/FPP/micefig2/4dpi.t.subset.marker.dotplot.pdf',
            width = 60)

publish_pdf('mission/FPP/micefig2/E.faecalis.4dpi.T.marker.pdf', width = 60)

## violin --------
sobj.tc <- sobj.tc |>
  filter(str_detect(orig.ident, 'SD')) |>
  mutate(group = str_remove(orig.ident, '-SD') |> fct_relevel('PBS'))

sobj.tc |> VlnPlot('Ifng', group.by = 'group')  
  
gh17 <- sobj.tc |>
  filter(type.fine == 'gd Th17 cells') |>
  bill.violin(c('Il17a','Il17f'), group.by = group, facet.ncol = 1) +
  labs(x = 'Group', fill = 'Group', y = 'Expression level',
       title = 'gd Th17 cells')

gh2 <- sobj.tc |>
  filter(type.fine == 'Th2 cells') |>
  bill.violin(c('Il4','Il6'), group.by = group, facet.ncol = 1) +
  labs(x = 'Group', fill = 'Group', y = 'Expression level',
       title = 'Th2 cells')

gcd8 <- sobj.tc |>
  filter(type.fine %in% c('CD8 T cells')) |>
  bill.violin(c('Ifng','Gzmb'), group.by = group, facet.ncol = 1) +
  labs(x = 'Group', fill = 'Group', y = 'Expression level',
       title = 'CD8 T cells')

gh2 + gh17 + gcd8 + plot_layout(guides = 'collect') &
  scale_fill_hue(direction = -1) &
  theme_jpub

publish_pdf('mission/FPP/micefig2/4dpi.t.subset.marker.violin.pdf', width = 70)

gh2$data |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  select(-.cell) |>
  write_csv('th2.4dpi.violin.csv')

gh17$data |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  select(-.cell) |>
  write_csv('th17.4dpi.violin.csv')

gcd8$data |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  select(-.cell) |>
  write_csv('cd8t.4dpi.violin.csv')

## save rds -----------
sobj.tc |> write_rds('mission/FPP/zhang21sa/zhang21tcell.rds')

sobj.tc <- read_rds('mission/FPP/zhang21sa/zhang21tcell.rds')
